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Abstract 

We report the results of a study of multiphase flow in porous media. A Darcy's law for steady 
multiphase flow was investigated for both binary and ternary amphiphilic flow. Linear flux-forcing 
relationships satisfying Onsager reciprocity were shown to be a good approximation of the simulation 
data. The dependence of the relative permeability coefficients on water saturation was investigated 
and showed good qualitative agreement with experimental data. Non-steady state invasion flows were 
investigated, with particular interest in the asymptotic residual oil saturation. The addition of surfactant 
to the invasive fluid was shown to signiflcantly reduce the residual oil saturation. 

1 Introduction 

The study of fluid flow in porous media is of great industrial importance. Applications include enhanced oil 
recovery, aquifer purification, containment of toxic and nuclear waste, geological flows of magma, chemical 
reactions in catalysts, and the study of blood flow through capillaries. The fluids of interest are frequently 
'multi-functional', for example, certain oil field fluids may be required not only to displace oil and gas and 
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transport them to the surface, but also to act as coolants and lubricants for the drill bit and for the transport 
of cuttings from within the well. Blood transports salts and sugars in solution, as well as white and red blood 
cells in colloidal suspension. Such fluids clearly require a more sophisticated description than that provided 
by single phase continuum approaches. Subjects of theoretical interest include the relationship between 
the macroscopic transport coefficients of porous media and their microscopic geometry , the geometrical 
properties of interfaces and the relationship between the transport coefficients and the morphology of 
the flow §,0. 

The flow of multiphase fluids in arbitrary geometries remains a challenge for conventional modelling 
techniques ^. In a top-down computational fluid dynamics (CFD) approach one would use a finite- 
difference or similar method to numerically integrate the Navier-Stokes equations with appropriate boundary 
conditions. Such methods become enormously complex when one wishes to simulate flow in an arbitrary 
geometry, or when one wishes to introduce the complex interfaces present in multiphase flow. The alternative 
atomistic approaches are not viable, as a 'typical' rock pore size is on the micron scale, far larger than the 
nanometre domain accessible with present day molecular dynamics. 

In order to make progress in the study of fluid flow in porous media, it is therefore necessary to simplify 
either the rock geometry or the fluid flow. Reducing the dimensionality of the flow enables larger systems 
to be studied computationally. The results from these simulations can then be compared with experimental 
results from micromodels pO[ |. In three dimensions, one can construct a simplified model of a single pore, 
and model the medium as a network of such simplified pores. The flow of a multiphase fluid through such 
a network is modelled based on some simplifying assumptions about the flow in each pore. Such models 
have been used extensively to study steady state flows [ pl| . A variation on the network model, 'invasion 
percolation', has been used to investigate the development of interfacial fronts in the invasive flows of interest 
in oil field situations ■ 

In recent years X-ray micro-tomography techniques have developed sufficient resolution to provide accu- 
rate three-dimensional digitizations of real rock structures. It is evidently preferable to study flow processes 
taking place in such realistic porous media. In order to do this a simplifled, computationally efficient model 
of multiphase fiow is required. Mesoscalc fluid models present clear advantages over conventional CFD 
techniques for this application. 

Lattice-gas and lattice-Boltzmann models have been used to study flow through porous media in both two 
and three dimensions. In two dimensions, lattice-gas models have been used to study both binary immiscible 
and ternary amphiphilic flow |]l^-[l5|l . In three dimensions the Rothman-Keller model has been used to study 
binary immiscible flow through a digitized slab of Fontainebleau sandstone ||^ . Studies of binary immiscible 
flow in three dimensions were also performed using a lattice-Boltzmann model [|6| |l7| . 

In this paper the three dimensional amphiphilic lattice gas model, described in detail in JTsf is applied 
to the study of flow through porous media. The adaptations of the model to handle arbitrarily complicated 
boundaries and fluid forcing are described. An investigation of steady state flows and comparison with 
Darcy's law is then presented for both two-phase binary immiscible and three-phase ternary amphiphilic 
flow. The simulation results are discussed in the context of linear nonequilibrium thermodynamics and 
recent extensions of Onsager theory to multiphase flow in porous media [p^-p3[. The behaviour of non- 
steady-state, 'invasion' flow is then simulated. Such simulations are intended to reproduce oil fleld flows, 
and the effect of surfactant on flow morphology and oil extraction is investigated. This work represents an 
extension of previous results obtained with the two-dimensional lattice gas model 
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2 The lattice gas model 



Our latticc-gas model is based on a microscopic, bottom-up approach, where dipolar amphipile particles are 
included alongside the immiscible oil and water species. Lattice-gas particles can have velocities c^, where 
1 < i < 6, and b is the number of velocities per site. We shall measure discrete time in units of one lattice 
timestep, so that a particle emerging from a collision at site x and time t with velocity will advect to site 
X + Ci where it may undergo the next collision. We let nf (x, t) e {0, 1} denote the presence (1) or absence 
(0) of a particle of species a G {R,B,A} {R, B, A denoting red (oil), blue (water) and green (amphiphile) 
species respectively) with velocity c^, at lattice site x G £ and time step t. The collection of all nf (x, t) for 
1 < i < b will be called the population state of the site; it is denoted by 

n(x,i)eA/' (1) 

where we have introduced the notation TV for the (finite) set of all distinct population states. The amphiphile 
particles also have an orientation denoted by (Ti (x, t). This orientation vector, which has fixed magnitude tr, 
specifies the orientation of the amphiphile particle at site x and time step t with velocity c^. The collection 
of the b vectors cr^ (x, t) at a given site x and time step t is called the orientation state. We also introduce 
the colour charge associated with a given site, 

g,(x,t) = nf (x,t)-nf (x,i), (2) 

as well as the total colour charge at a site, 

b 

Q (x, t)=^q, (x, t) . (3) 

i=l 

The state of the model at site x and time t is completely specified by the population state and orientation 
state of all the sites. The time evolution of the system is an alternation between an advection or propagation 
step and a collision step. In the first of these, the particles move in the direction of their velocity vectors to 
new lattice sites. This is described mathematically by the replacements 

<(x + c„t + l) ^<(x,t), (4) 

o-i (x + Cj,t + 1) ^ (Tj (x, i) , (5) 

for all X G £, 1 < i < 6 and a G {R, B, A}. That is, particles with velocity simply move from point x 
to point X + Ci in one time step. In the collision step, the newly arrived particles interact, resulting in new 
momenta and surfactant orientations. The collisional change in the state at a lattice site x is required to 
conserve the mass of each species present 

b 

p"(x,t)^^<(x,0, (6) 
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as well as the D-dimensional momentum vector 



p(x,i) EE^^C,7ZnX,t), (7) 

a i 

(where we have assumed for simplicity that the particles all carry unit mass). Thus, the set TV of population 
states at each site is partitioned into equivalence classes of population states having the same values of these 
conserved quantities. 

We assume that the characteristic time for collisional and orientational relaxation is sufficiently fast in 
comparison to that of the propagation that we can model this probability density as the Gibbsian equilibrium 
corresponding to a local, sitewise Hamiltonian function; that is 

p(,')^lcxp[-/3Jf(,s')], (8) 

where (3 is an inverse temperature, H{s') is the energy associated with collision outcome s', and Z is the 
equivalence-class partition function. The sitewise Hamiltonian function for our model has been previously 
derived and described in detail for the two-dimensional version of the model p^, and we use the same one 
here. It is 

H{s') = J ■ (aE + ^P) + cr' ■ (eE + CP) + J ■ {(-£ + (V) + ^v(x, i)', (9) 
where we have introduced the colour flux vector of an outgoing state, 

b 

J(x,t) =^c,q,'(x,t), (10) 

1=1 

the total director of a site, 

b 

cr (x, i) = ^ cr, (x, i) . (11) 

4=1 

the dipolar flux tensor of an outgoing state, 

b 

^(x,t) = ^c,^Tax,t), (12) 



the colour field vector, 



the dipolar field vector. 



the colour field gradient tensor. 



E(x,t) =^c,g(x + c,;,i), (13) 



P(x,t) = -^c,^(x + c,;,t), (14) 



f (X,t) EEE^C,E(x-t-C„i), (15) 



4=1 
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the dipolar field gradient tensor, 



b 



(16) 



defined in terms of the scalar director field, 
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5(x,t) 



(17) 



1=1 



and the kinetic energy of the particles at a site. 




(18) 



where v is the average velocity of all particles at a site, the mass of the particles is taken as unity, and a, ^, 
e, C, and 5 are coupling constants. To maintain consistency with previous work we use the coupling constants 
as previously defined in ||l^ . The values of these constants are: 



These values were chosen in order to maximise the desired behaviour of sending surfactant to oil-water 
interfaces while retaining the necessary micellar binary water-surfactant phases. It should be noted that the 
inverse temperature-like parameter (5 (whose numerical value is varied in this paper) is not related in the 
conventional way to the kinetic energy. For a discussion of the introduction of this parameter into lattice 
gases we refer the reader to the original work by Chen, Chen, Doolen and Lee and Chan and Liang [p6| . 
Eqs. d) -(|l6|) were derived by assuming that there is an interaction potential between colour charges, and 
that the surfactant particles arc like "colour dipolcs" in this context ||2^. The term paramctcriscd by a 
models the interaction of colour charges with surrounding colour charges as in the original Rothman-Keller 
model |27t ; that parameterised by ^ describes the interaction of colour charges with surrounding colour 
dipoles; that parameterised by e accounts for the interaction of colour dipoles with surrounding colour 
charges (alignment of surfactant molecules across oil- water interfaces); and finally that parameterised by 
Q describes the interaction of colour dipolcs with surrounding colour dipoles (corresponding to interfacial 
bending energy or "stiffness"). This model has been extensively studied in two dimensions p^, p8|-^, 
and the three-dimensional implementation employed in the present paper is described in more detail by 
Boghosian, Coveney and Love jl8|. 

2.1 Porous media 

Porous media are introduced into lattice-gas simulations by adding another bit to the description of the 
state at a lattice site. This bit takes the value one at rock sites and zero at pore states. In order to 
realistically simulate oil field multiphase flows it is necessary to include the wetting properties of the porous 
media. Within the lattice-gas model described above, the rock sites are assigned a colour charge w which 
may range from —26 < w < 26. Thus w = —26 represents a rock which is completely water wetting. 
This implementation of rock sites is identical to that used in the two-dimensional implementation of the 



The porous medium used for all ensuing simulations is a subset of a Fontaineblcau sandstone X-ray 
microtomography dataset. The original dataset was obtained by Exxon Mobil research using the synchrotron 



a = 1.0, e = 2.0, = 0.75, C = 0.5 



model |T|,[5). 
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source at Brookhaven. The total dataset consisted of approximately 19.5 million voxels at a scale of 7.5^m, 
over a volume roughly 2mm on a side. This dataset was also used by Olson and Rothman j^]. The porosity 
(ratio of void space to total volume) of the medium is 0.58. Fontaincblcau sandstones are widely used 
as they are one of the cleanest sandstones available, consisting mainly of reasonably monodisperse quartz 
sand grains. They also contain low concentrations of paramagnetic iron salts, minimising the difficulties 
associated with flow imaging by techniques such as nuclear magnetic resonance. The system size chosen for 
all the following simulations was 64'^. This is not the largest system size attainable, but was chosen as a 
reasonable size to allow large numbers of simulations to be performed as we are interested in the behaviour 
of our model over a wide range of fluid composition and fluid forcing levels, as well as requiring considerable 
ensemble averaging to obtain reasonable precision on our measurements. 

2.2 Fluid forcing 

Fluid forcing is implemented in the code in order to reproduce the effect of one fluid experiencing a buoyancy 
effect due to gravity. Forcing is carried out by selecting random sites across the lattice and adding or removing 
momentum as required. Momentum is added until a fixed amount of momentum per lattice site of the forced 
fluid has been added. The forcing level referred to below is the average amount of momentum added per 
forced fluid site. Forcing in is in the positive z-direction. 

2.3 Fluid-solid boundary conditions 

Two types of boundary conditions can be used at fluid-solid boundaries: slip and no-slip conditions. In the 
slip condition, a particle incident on the solid wall undergoes specular reflection, that is, the particle preserves 
its component of momentum parallel to the wall, and has its component of momentum perpendicular to the 
wall reversed. In no-slip boundary conditions, particles incident on the wall have their momentum reversed, 
leading to the alternative description of 'bounce-back' conditions. In CFD simulations the no-slip condition 
is typically used, in that case by imposing a zero tangential fluid velocity at the boundary. 

No-slip boundary conditions may be trivially implemented for lattice-gas and lattice-Boltzmann simu- 
lations. Particles (or single particle distribution functions in the case of lattice-Boltzmann) which advect 
onto a site occupied by a rock particle are moved into states with the opposite momentum. This will cor- 
respond to a zero velocity condition between the lattice sites adjacent to the boundary and the boundary 
itself However, if the fluid has a wetting interaction with the solid, such a boundary condition may 
not be appropriate. Recent non-equilibrium molecular dynamics simulations have illustrated this, showing 
that a significant slip velocity may exist for such situations pBf . Mixed boundary conditions allowing for a 
non-zero slip velocity at the wall could be introduced into the lattice-gas model by including more complex 
fluid-solid collision rules. However, this possibility remains unexplored, and we do not pursue it in the 
following simulations. 

2.4 Simulation box boundary conditions 

A typical porous medium is not periodic, and so periodic boundary conditions at the edge of the simulation 
box cannot be used. It is important to keep some periodicity in the direction of the flow to maintain 
continuity of the flow of fluids. In order to achieve periodicity in the direction of the flow, a mirror image 
of the porous medium is built in the direction of the flow. The use of the mirror image ensures that all the 
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pores are connected, and so periodic boundary conditions can be used in this direction. This method was 
first used by Olson and Rothman in their three dimensional lattice gas studies 0,0- For the boundary 
conditions in the directions transverse to the flow it is not practical to use the same method because it 
increases the box size by a factor of two for each direction. 

There are two possible boundary conditions that can be applied: 'wall conditions' or 'buffer conditions'. 
The wall conditions consist of creating a plane of obstacle sites at the simulation box edge in the transverse 
directions. Every particle coming to the wall is normally bounced back. The buffer conditions consist of 
removing a region of width two lattice sites around the faces of the box and then placing a wall around 
the faces. Particles can then flow in and out of this reservoir into the pore space. The different boundary 
conditions were tested on a simulation of a one component ffuid ffowing in a 32'^ porous medium. A reference 
value for the permeability has been chosen using a 64^ porous medium and calculating the flow only in the 
core (assuming that the boundary conditions do not influence the flow in the core of the material) . Using the 
wall conditions the permeability is half the reference permeability, while using the buffer conditions, a value 
close to the reference value was found. The buffer conditions have therefore been used for all subsequent 
simulations. 

For invasion simulations the boundary conditions in the flow direction must simulate the effect of the 
system being connected to an inflnite reservoir of invasive fluid at the 'bottom' of the system. At the 'top' 
of the system the boundary conditions must reproduce the effect of a reservoir of fluid displaced from the 
porous media. For the binary case where the invasive ffuid is water penetrating a medium saturated by oil 
the boundary conditions arc relatively simple. Oil particles leaving the top of the system are transferred to 
the bottom and recoloured, so that an equal number of water particles enter the bottom of the system. The 
total number of particles in the simulation is conserved, although the composition of the system obviously 
changes. 

For the ternary invasion case the situation is more complex. The boundary conditions described above 
create an artificial oil-water interface at the system boundaries in the flow direction. The surfactant particles 
show a strong tendency to adsorb at these 'interfaces'. This gives rise to two simulation artefacts. Firstly, 
adding a constant amount of surfactant to the system at each timestep does not give a water: surfactant ratio 
in the invasive fluid equal to that added. This is because the surfactant particles adsorb to the fictitious 
oil- water interface at the bottom of the system rather than entering the porous medium. The second spurious 
effect occurs when the invasive fluid front reaches the top of the system. The 'oil reservoir' simulated by 
the upper boundary condition is then in contact with a high concentration water-surfactant solution. This 
water surfactant solution then solubiliscs oil out of the oil reservoir above the system, leading to an increase 
in the oil concentration in the system. Clearly both these effects are highly undesirable. 

The first artefact can be removed by actively controlling the amount of surfactant injected to compensate 
for the loss of surfactant to the 'interface' at the bottom of the system. The second artefact is also simply 
removed: one replaces the oil reservoir by a reservoir with the same composition as the top layer of sites in 
the system. 
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3 Steady-state flow - Darcy's law 



We validate an extension of the single phase Darcy's law to multiphase flow which explicitly admits a viscous 
coupling between fluid phases. The single phase Darcy's law relationship is: 

J = --(Vp-pg), (19) 
M 

where J is the fluid flux, k is the permeability, /i is the kinematic viscosity, and (Vp — pg) is the force on 
the fluid due to pressure gradients and gravity. 

As previously noted there is no widely accepted hydrodynamic description of multiphase fluids. Similarly 
there is no universally agreed extension of Darcy's law for multiphase flow. The single phase Darcy relation 
has been commonly extended as: 

J, = ^A:y(5)-X„ (20) 

where fcy (5*) is the relative permeability coefRcient (depending only on the saturation), Jj is the flux of the 
ith species and Xi is the body force acting on the jth component. 

A number of questions arise concerning these phenomenological descriptions. Firstly, what is the domain 
of validity of the linear relationship between flux and forcing? As the flow is still governed by the (non- 
linear) Navier-Stokes equations the linear dependence in (^9|) and ( pO[ ) must break down at some forcing 
value. Secondly, the dependence of the permeability coefficients on the rock microgeometry is of considerable 
interest. Ideally one would like to determine the coefficients k and kij independently, and then use eqns ( [l9| ) 
and ( ^o|) to predict flow rates. Finally any symmetry properties of the matrix kij are of considerable interest 
from a theoretical standpoint. 

Equation (|20| ) has the form of a linear phenomenological law. Such descriptions of transport processes 
form the theory of linear nonequilibrium thermodynamics. The general theory of such processes was first 
described by Onsager in 1931 ||l9|,^. On the basis of the reversibility of the microscopic dynamics of atomic 
motion Onsager derived the reciprocity relations Q 

ki2 = fcai. (21) 



The validity of these relations additionally requires the thermodynamic forces Xj and fluxes J,; to be conju- 
gate in the sense that the entropy production can be written ■ X^. 

Flekk0y and Pride recently demonstrated that the fluxes and forces in Darcy's law are indeed re- 
ciprocal. By assuming a particular form for the mechanism of forcing they were able to demonstrate the 
validity of reciprocity in the infinite capillary number limit, where the interface between the fluids does not 
move. They went on to consider two additional cases, firstly where the interface relaxes to its equilibrium 
position from a small displacement, and secondly where surfactant is present in the system. The additional 
entropy production considered in the first case arises from the work done by the interface on the bulk fluid. 
In the second case surfactant gradients produce additional tangential forces on the interfaces. In these two 
more complex cases the fluxes and forces identified are somewhat sophisticated: for example the additional 

^ It is interesting to note that in his papers Onsager was unwiUing to commit himself as to whether the dynamics of atomic 
motion was actually irreversible. Writing only six years after Heisenberg's first paper on matrix mechnics, quantum mechanics 
was insufficiently well established as a reversible theory of atomic dynamics to settle this issue. 
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fluxes involved are the Fourier components of the displacement of the interface and surfactant concentration 
from their equilibrium values. Reciprocity was shown to hold in both cases. 

The domain of validity of these results and their applicability to results obtained from our simulations is 
clearly of great interest. Flekk0y and Pride point out that any local nonlincarity (such as bubble breakup) in 
the flow will invalidate Although the presence of such local nonlinear phenomena may result in an average 
flow behaviour which obeys a force-flux relationship of the form (|2^), there is no reason that the coefficients 
kij should be the same as those for completely linear hydrodynamics. The presence of reciprocity in the 
results obtained by Olson and Rothman Q and in the two dimensional implementation of our model fl^ 
to within the simulation accuracy is therefore still a complete mystery from the point of view of linear 
noncquilibrium thermodynamics. 

3.1 Single phase and binary immiscible flow 

We performed simulations which address the issues discussed above in the case of binary immiscible flow. 
Firstly we measured the flux as a function of forcing to investigate the validity of linear relationships such 
as (|l^) and (|2^). We then simulated single phase flow in order to obtain the reference permeability of our 
medium. Five independent simulations were performed for 2000 time steps on 16 processors, for forcing 
levels 0.005, 0.03, 0.05, 0.1, 0.15. The flux was time-averaged after allowing 500 timcsteps for the flow to 
reach a steady state. The results of these simulations are displayed in Figure |l|. 

The behaviour of the flux as a function of forcing was then investigated for binary immiscible flow at the 
same values of forcing as in the single phase case. Two types of simulations were performed: one in which oil 
was forced, and one in which water was forced. Two independent simulations were performed for 2000 time 
steps at each forcing level and for each forced fluid. The flux was time-averaged after allowing 500 timesteps 
for the flow to reach a steady state. The results of these simulations are displayed in Figure |l|. 




0.05 0.1 0.15 

Forcing IgvgI (lattice units) 



0.05 0.1 0.15 

Forcing level (lattice units) 



(a) 



(b) 



Figure 1 : Veriflcation of Darcy's law for a) Single phase and binary fluids when forced, b) Flux of binary 
fluids when unforced. Flux is normalised to show average momentum at each lattice site in the pore space. 
Diamonds show singl;e phase flux, squares show oil flux and circles show water flux. 



The squares correspond to the flux of oil when it is forced and the circles correspond to the flux of 
water. A linear variation of the flux versus the force is observed, over the whole range of forcing levels. The 
relative permeability coefficient is smaller in the case of water because of the wettability of the rock. One 
expects the water phase to flow adjacent to the rock and therefore experience a greater drag force than the 
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oil phase. Visualisation of the flow morphology for our simulations verifies this expectation. As can be seen 
from figure |l(b)| the reciprocity relation kij = kji is valid within the statistical uncertainty. 

These results are consistent with those obtained from the two dimensional implementation of our model [ p^ 
The nonlinear deviations from Darcy's law observed at high forcing in two dimensions was not observed in 
three dimensions. The magnitude of the cross terms in three dimensions is reduced by an order of mag- 
nitude. This is consistent with previous three dimensional results, and the magnitude of the cross terms 
in [p^ was believed to be due to the dimensionality of the model. In two dimensions our model was shown 
to exhibit a capillary threshold at low forcing. The flux of the unforced fluid would only become non-zero 
at some threshold forcing level. Such effects are not observed in the results presented above, however, the 
porosity of the medium used here (0.58) is close to that of the medium with the largest porosity used in 
two dimensions (0.60). The capillary threshold in two dimensions systematically decreased with increasing 
porosity. The considerable extra computational cost of performing calculations in three dimensions did not 
permit a systematic study of the influence of the geometry of the porous medium on Darcy's law behaviour. 

The behaviour of the matrix kij as a function of the fluid composition is of considerable interest, and 
enables us to connect the lattice-gas model to previous experiments and network models in which this 
dependence has been investigated 3^. In order to perform simulations over a range of compositions 
within available computational resources the forcing level was kept constant at 0.05. Simulations were 
performed at reduced density of water 0.05 < Pu, < 0.45 at intervals of 0.05. The total reduced density was 
held constant at 0.5. Two types simulations were performed, one in which oil was forced, and one in which 
water was forced. Two independent simulations were performed for 2000 time steps at each composition and 
for each forced fluid. The flux was time averaged after allowing 500 timestcps for the flow to reach a steady 
state. The results of these simulations are displayed in flgure 2(a) and ^ 

Figure || shows the diagonal components of kij for the experiments conducted by Wyckoff and Botset . 
This data was obtained from flows of water and carbon dioxide in unconsolidated sand packs. The data 
displayed is for the sand specimen with porosity 0.47, the closest porosity to that of our virtual Fontainebleau 
sandstone. 

The general features of these graphs can be illuminated considerably by visualisation of the flow mor- 
phology in our simulations. For low water (wetting fluid) saturations small bubbles of water exist in a 
background of oil. These bubbles are predominantly adjacent to the rock as one would expect. Increasing 
the water saturation increases the coverage of the rock by water, and pores filled with water appear. The 
morphology of the oil phase is fully connected through the pores until the water saturation reaches 0.8. For 
saturation greater than 0.8, there is no connected path of oil through the porous medium. 

How do the changes in oil morphology relate to the features seen in figures || and ^ First consider the 
diagonal components of kij. For low water saturations the flux of water is zero below some critical saturation, 
around 0.2 for our simulation data and 0.4 for the comparison experimental data of | |35[ |. This can be ascribed 
to the presence of many small droplets of wetting fluid being adsorbed onto the media and therefore unable 
to flow. The larger threshold in the experimental case can be interpreted as a consequence of the differences 
in forcing between the two cases. The water was driven in the experiments by a pressure gradient, and so 
the net momentum transferred to the fluid will depend on the morphology. In particular, the fluid will not 
flow until a connected path exists across the medium. This is not the case in our simulations where the 
forcing reproduces the effects of the wetting fluid being driven by buoyancy forces. 

A similar argument may be applied to the reduction of the non-wetting fluid flux to zero at a saturation 
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0.2 0.4 0.6 0.8 1 

Water saturation (lattice units) 



(a) 

Figure 2: Variation of the diagonal components of the relative permeability matrix fcy with wetting fluid 
saturation. Lattice-gas simulation data for a 64'^ system using the porous medium described above (lines). 
The medium has porosity 0.58. Experimental data on C02/water flow through unconsolidated sand sample 
of porosity 0.47 from Wyckoff and Botset (symbols) |^^. Squares show gas permeability (non-wetting fluid), 
circles show water permeability (wetting fluid). 




0.4 0.6 
Water saturation 



Figure 3: Variation of the off-diagonal components of the relative permeability matrix kij with wetting fluid 
saturation. The squares show relative permeability for oil (non-wetting fluid) when water is forced, the 
circles show relative permeability for water (wetting fluid) when oil is forced. 
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of around 0.9. At this saturation the non-wetting fluid ceases to be connected. It should be pointed out that 
both these effects are illustrative of a capillary threshold for the flow. Such a threshold was not observed in 
our study of the behaviour of flux versus forcing above. However, for the composition used in that study 
both fluids were continuously connected across the full dimensions of the porous medium. The agreement 
between the permeability of water in figures |^ and |^ is quite good. However, below a water saturation 
of 0.3 the non-wetting permeabilities are quite different. This presumably arises because the experimental 
data is obtained for water and carbon dioxide, i.e. a liquid/gas flow where the fluids have very different 
hydrodynamic behaviour. Our immiscible fluids are hydrodynamically identical. An extensive study of the 
dependence of fcy on the viscosity ratio of the two fluids was performed in |pl| , and features such as that 
reported in figure ^ below a water saturation of 0.3 were obtained. 

The statistical uncertainties on the off-diagonal components of kij make it impossible to formulate any 
definitive statements; nevertheless a few comments are appropriate. Firstly we see that, within the statistical 
uncertainties kij is a symmetric matrix. As pointed out above this is not a consequence of Onsager reciprocity, 
and future work should provide a more rigorous test of the symmetry properties of kij in the non-linear flow 
regimes accessible to our model. Secondly, the maximum flux of the unforced fluid occurs at a water 
saturation of 0.8, that is, at the percolation transition for the unforced fluid. This is consistent with the 
intuitive picture that viscous coupling is reduced for a flow of bubbles in a background phase. 

These results show the same qualitative connection between connectivity of the oil and water phases and 
the features of the kij dependence on saturation. The signiflcant exception is the absence of lubrication 
effects in three dimensions. The oil (non-wetting) flux at low water saturation was significantly larger than 
that of a pure oil phase in two dimensions [Q . No such effect was observed here. 

3.2 Ternary amphiphilic flow 

In this section the infiucncc of surfactant on the linear flux-force description is investigated. We retain the 
form of the multiphase Darcy's law presented above. In this case the relative permeability matrix kij is 
a 3 X 3 matrix which explicitly allows coupling between any pair of species. The influence of surfactant 
will be described with reference to the questions addressed above, namely a) When is the linear flux-force 
relationship valid for ternary flow? b) What symmetry properties can be demonstrated for the matrix fcy ? 
c) What is the influence of surfactant on the magnitudes of the diagonal and off-diagonal components of 
kij as compared with the binary case? In the two dimensional version of the model the significant effect 
of surfactant was to reduce the capillary threshold by reducing the oil-water intcrfacial tension. We do not 
observe a capillary threshold in the binary Darcy's law simulations and so perform simulations to investigate 
the influence of surfactant on the general features of the binary Darcy law behaviour discussed above. 

We performed simulations for 2000 time steps for 5 forcing values. Two types of simulations were 
performed, one in which oil was forced, and one in which water was forced. Four values of water reduced 
density were investigated: 0.05, 0.15, 0.35 and 0.45. The surfactant reduced density and total reduced 
density were kept constant at 0.1 and 0.6, respectively. Fluxes were averaged over for 1500 time steps after 
allowing 500 timesteps for the flow to reach a steady state. Because simulations in which the surfactant is 
forced were not performed the ki3,k23 and fcsa components of kij were not investigated in these simulations. 
The kii and ^22 components of kij are plotted as a function of water saturation in figure |[ These diagonal 
components are considerably reduced as compared with the binary case. 

The infiuence of the surfactant on the fluid-fluid coupling may be investigated by examining the behaviour 
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Figure 4: Variation of the kn (squares) and ^22 (circles) components of the relative permeability matrix ki 
with water saturation for ternary amphiphilic steady-state flow. 
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of the off-diagonal components of the relative permeability matrix. For low water saturations the measured 
fluxes of the unforced fluid are swamped by the statistical errors. A comparison of the flux-forcing relationship 
for water saturation 0.7 is shown in figure 5(a)| . Within the statistical errors this sector of the relative 
permeability matrix appears to remain symmetric in the presence of surfactant. In figure [5(b)] the surfactant 
flux when oil or water are forced is shown, as a function of forcing. This validates the extension of Darcy's 
linear relationship, eqn ( pO| ) to the three phase case, at least for ternary amphiphilic flows. A linear response 
of the surfactant to oil forcing is also seen for water saturation 0.05 and 0.15. 



Figure |5(b)| implies that, within statistical uncertainty, fcsi = ^32 at this water saturation. This deserves 
some consideration. Clearly, any postulated symmetry of fcy does not imply k^i = k^2- The microdynamics 
of our model is invariant under interchange of oil and water particles. However, because the wettability of 
the rock as defined above is —26, macroscopic parameters such as fc^ are not symmetric with respect to 
interchange of oil and water. If most of the surfactant is adsorbed to interfaces in the flow far from solid-fluid 
boundaries, one might expect the water-surfactant and oil-surfactant coupling to be symmetric. Without 



accurate data for a range of water saturations it is of course possible that the agreement displayed in p(b) 
is coincidental. 




0.05 0.1 0.15 

Forcing level (lattice units) 




0.05 0.1 0.15 

Forcing level (lattice units) 



(a) 



(b) 



Figure 5: Fluid-fluid coupling for steady state ternary amphiphilic flow. System size 64"^. a) Oil flux (squares) 
and water flux (circles) when water and oil are forced, respectively. Fluxes measured for a water saturation 
of 0.7. b) Surfactant flux when oil is forced (squares), and when water is forced (circles). 



Visualisation of the flow morphology for these simulations confirms that the surfactant is behaving as 
expected. The preference of surfactant to adsorb at oil-water interfaces leads to an additional effect for 
the water wetting medium studied here. An oil-solid boundary plays the same role as a water-surfactant 
interface, and so we expect surfactants to adsorb to the rock in such cases, thereby modifying the wetting 
properties of the rock. Such behaviour has been observed in our simulations. The surfactant can take 
on three roles: firstly it can exist in bulk oil or bulk water as (at this concentration) wormlike micelles. 
When the saturation of one fiuid is much greater than the other this is the predominant behaviour. The 
micelles are carried along in the fiow, leading to the water-surfactant and oil-surfactant coupling discussed 
above. Secondly the surfactants adsorb to oil-water interfaces as expected. Thirdly the surfactant adsorbs 
to the porous media. For intermediate saturations all the oil flows either in pores coated by a thin layer 
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of surfactant, or if water coats the rock surfaces a layer of surfactant is adsorbed to the oil- water interface. 
The tendency of the surfactant to enable oil to flow adjacent to the rock, rather than being lubricated by a 
layer of water, may be one factor contributing to the reduction of the diagonal components of fcy discussed 
above. 



4 Imbibition Simulations 

The steady state studies in the previous section attempt to characterise the flow in terms of the simple Darcy's 
law phenomenological description. Such phenomenological relationships are frequently used as constitutive 
relationships in coarse-grained continuum models. These models are used to study the non-steady state 
flows of interest in oil-rcscvoir flows and ground water pollution remediation. Such an approach could be 
described as ad hoc, at best. One advantage of the lattice-gas technique is that the model provides a general 
description of the flow, valid for both steady and non-steady flows. 

Potentially one of the most commercially rewarding applications of surfactants is their use in enhanced 
oil recovery. In oil flcld flows, oil is displaced and carried to the surface by water. Two limiting factors exist 
for this technique. Firstly, the ratio of oil to water must be high enough for the process to be commercially 
viable. Secondly, the total amount of oil cxtractablc is limited. As one might expect, the immiscibility 
of oil and water eventually confounds efforts to solubilise one into the other. These limitations mean that 
most commercially developed oil fields have between 40% and 60% of the available oil remaining after 
primary and secondary recovery operations. Addition of surfactant to the extraction fluid should increase 
the solubilisation of oil into the water, thereby increasing well productivity. It should be noted that the 
amounts of surfactant needed and their associated cost mean such techniques will not flnd application until 
the world oil price has increased substantially. 

We perform simulations intended to reproduce such oil fleld extraction flows. The porous media is initially 
saturated with oil at a reduced density of 0.5, and a water surfactant mixture is forced into one end of the 
media. The subsequent flow morphology and time history of oil concentration remaining in the media are 
then studied. The wettability of the porous media is kept constant —26, i.e. completely water wetting. 
Invasion into such a porous medium is referred to as imbibition, from the medium's natural tendency to 
imbibe the wetting fluid even in the absence of forcing. These simulations extend work performed in two 



dimensions by Fowler and Coveney and Maillet and Coveney |14, [5[. 

We vary two parameters in the invasion study, the fluid forcing level and the fraction of surfactant in the 
invasive fluid. The range of input forcing is limited by two constraints. The forcing must be low enough that 
the underlying lattice-gas dynamics still reproduces Navicr-Stokes behaviour, and must be high enough that 
the simulation reaches the percolation point within simulation lengths attainable with current resources. 

We first study the extraction of oil by an invasive fluid comprising pure water. Five independent simu- 
lations were performed for each forcing level at forcing levels 0.01, 0.02, 0.03, 0.04, and 0.05. The temporal 
evolution of the oil saturation is shown in flgure |[ There are three production regimes observed here. By 
visualising the densities of oil and water we may identify these as follows: an initial slow regime prior to 
water percolation, a regime subsequent to water percolation and prior to oil depercolation in which the 
establishment of high flow rate channels causes rapid oil production, and a third regime after the end of oil 
percolation in which the oil saturation approaches its asymptotic value. 

In the two dimensional version of our model the flrst and last of the extraction regimes described above 
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were observed. This may be a dimensional effect, as the second regime in which oil is extracted rapidly by 
water domains connected across the medium requires both oil and water domains which extend across the 
medium. Such a bicontinuous morphology of oil and water in the porous media can only occur in three 
dimensions. 

3e+06 I ■ 1 ■ 1 ' 1 ■ 1 ■ 1 




I ■ ' ■ ' ' ' ■ ' ■ 1 

2000 4000 6000 8000 10000 

Time (time steps) 

Figure 6: Number of oil particles in system as a function of time for various forcing levels in a binary system. 
The forcing level is from right to left 0.01, 0.02, 0.03, 0.04, 0.05 (lattice units). Error bars show the standard 
deviation on the mean of five independent simulations. 

As pointed out by Maillot and Covcney, the third regime here is of most interest in the context of enhanced 
oil extraction. In order to investigate the influence of surfactant on the flow morphology and prc-pcrcolation 
behaviour 5 independent simulations were performed for each value of forcing and invasive composition, for 
the 5 forcing levels used in the binary case, and for invasive fluid water:surfactant ratios of 4 : 1 and 2:1. 
These are ratios which lead to equilibrium phases consisting of wormlike micelles, droplet phases and sponge 
microemulsions respectively. It is expected that the flow conditions and porous media will cause significant 
deviations from the equilibrium morphology. The time evolution of oil saturation as a function of surfactant 
concentration is shown in figure ^. 

This graph has a number of notable features. Firstly, the production rate prior to water percolation is 
reduced by the presence of surfactant. We attribute this to the reduction of surface tension caused by the 
surfactant, which leads to a more convoluted interface, rather than the 'piston-like' displacement we observed 
in the binary invasion case. However, the surfactant does systematically reduce the time required for water 
percolation as a consequence of the more convoluted interface. The reduction in oil production relative to 
the binary system at early times in the ternary system was also reported for the two dimensional version of 
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Figure 7: Number of oil particles in system as a function of time for various invasive fluid compositions in a 
ternary system. Error bars show the standard deviation on the mean of five independent simulations. 
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our model by Coveney et al.. Only by performing simulations which reached longer time scales could Maillet 
and Coveney assess the impact of surfactant on the residual oil saturation . 

In order to assess the impact of surfactant on the asymptotic oil saturation further simulations were 
performed for forcings 0.01,0.02,0.03 and 0.04. Invasive fluid water-surfactant ratios of 4 : 1 and 2 : 1 and 
binary invasion without surfactant were simulated for 20000 timcstcps. The oil saturation reaches a steady 
state value after 15000 timesteps. The residual saturation was averaged over the last 5000 timesteps and a 
statistical error derived from the standard deviation on the mean. The results are plotted in figure g. 
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Figure 8: Residual oil saturation (%) for binary (circles), water:surfactant ratio 4 : 1 (squares) and 2 : 1 
(triangles). 



Figure ^ bears out our supposition that addition of surfactant to the extraction fluid should have a 
significant impact on the residual oil saturation in the media. In fact, the residual oil saturation is reduced 
by a factor of five between the binary case and the case with the highest water-surfactant ratio in the 
extraction fluid. This is illustrated dramatically by visualising the oil saturation in the medium at time step 
20000 for the three types of invasive fluid, as shown in figure 

These results are in complete agreement with those obtained in two dimensions by Maillet and Coveney ||l5|] 
In the two dimensional case the enhancement of oil recovery was shown to be due to the emulsification of 
the residual oil into the extraction fluid as surfactant coated droplets. In the next section we visualise 
the behaviour of the surfactant and compute structural characteristics of the interface to investigate the 
mechanisms involved in recovery enhancement in three dimensions. 



18 




Figure 9: Isosurfacc of residual oil at timestep 20000. a) No surfactant, b) Water-surfactant ratio 4:1. c) 
Water-surfactant ratio 2:1. Isosurface shows oil concentration at a value of 5 particles per site. Forcing 
level 0.01. Porous media is not displayed for clarity. 

5 Interfacial Morphology 

To better understand the role of surfactant in enhanced recovery we analyse the behaviour of the oil-water 
interface prior to invasive fluid percolation as a function of surfactant concentration. Visualisation of the 
interface shows that it is much more convoluted in the ternary case. Maillet and Coveney computed the 
fractal dimension of such interfaces in the two dimensional version of our model. This fractal dimension was 
largely determined by the fractal dimension of the porous medium itself. To avoid the additional complexity 
of such calculations in three dimensions we compute the interfacial distribution function: 

where /i(x) takes the value one if the site at x is at the interface and takes the value zero otherwise. Interfacial 
sites are defined as those at which the colour charge is zero. This function is useful as it allows us to define 
precisely the interfacial position: 

z = 5]/(z)z. (23) 

Prior to percolation the interface moves with constant speed through the system in both binary and 
ternary simulations. The interfacial speed, defined as z, is plotted in figure ^ The speed for low or inter- 
mediate surfactant concentration is proportional to the forcing level as expected. This result is equivalent to 
the linear dependence of early time extraction rate on forcing level observed in the two dimensional imple- 
mentation of our model. For high surfactant concentration the interfacial speed appears to be independent 
of forcing, a somewhat suprising result. From direct visualisation of the fluid interfaces it is clear that the 
picture of a well-defined interface moving coherently through the porous media is inaccurate for high sur- 
factant concentration. One anticipates that the surfactant will cause significant growth in the width of the 
interface as a result of emulsification processes. In figure ^ we show the interfacial distribution function 
as a function of position relative to z for binary fluid invasion and for water-surfactant ration of 1 : 2. All 
values of forcing are plotted. 

In the binary case significant growth of the interface does not occur over the timescalc of the simula- 
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tion. The addition of surfactant clearly causes the interface to grow significantly. The independence of 
the interfacial speed and the forcing level arises because emulsification processes dominate prior to water 
percolation. 
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Figure 10: Interfacial speed z for binary (circles), watcr:surfactant ratio 4 : 1 (squares) and 2 : 1 (diamonds). 



6 Conclusions 

The three dimensional lattice gas model for the non-equilibrium dynamics of amphiphilic fluids has been 
extended to simulate multiphase flow in porous media. The description of macroscopic behaviour of binary 
immiscible steady state flows by an extension of Darcy's law with an explicit viscous coupling between species 
has been verified. The extension of this Darcy's law description to ternary amphiphilic fluids has also been 
investigated 

In the binary case the dependence of the relative permeabilities for the forced fluids bear a strong similarity 
to results obtained from network models and flows in sands ||Tl|,Q. The cross coefficients satisfy a reciprocal 
relationship ki2 = fei, within the statistical uncertainty of the lattice-gas model. As pointed out by Olson 
and Rothman in 0, and this is not a consequence of Onsager reciprocity, but is nevertheless now a 

reasonably well established relation without theoretical foundation. The cross-coefficients are much smaller 
than those observed in the two-dimensional realisation of our model, as expected. 

In the ternary case the linear flux-forcing relationship has been verified for intermediate wetting fluid 
saturations. The relative permeabilities are reduced by the presence of surfactant. The presence of an 
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Figure 11: Scaled interfacial distributions for binary invasion (squares) and invasive water:surfactant ratio 
2 : 1 (points). 
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additional species increases the difficulty of studying these systems. It should be noted that a systematic 
study of Darcy's law behaviour as a function of surface tension is possible in the model by varying the inverse 
tcmpcrature-likc parameter (3. Such a study could elucidate which aspects of the ternary behaviour are due 
to the modification of surface tension by surfactant and which arc due to other factors, such as modification 
of the wetting properties of the media. 

Invasion simulations performed with both binary and ternary amphiphilic fluids show three regimes for 
production of the defending fluid (oil). Prior to invasive fluid percolation the oil is produced slowly, after 
percolation and prior to the end of oil percolation rapid production occurs, followed by a third regime in 
which the oil saturation approaches an asymptotic residual value. Addition of surfactant to the invasive fluid 
reduces the asymptotic residual oil saturation by a factor of five. The pre-percolation interfacial behaviour 
is also strongly affected by the presence of surfactant. The interfacial behaviour changes from a 'piston-like' 
movement of an interfacial region of fixed width to a motion dominated by the growth of the interfacial 
region. The contributions to the behaviour of capillary fingering and other surfactant-induced emulsification 
processes could be elucidated by a systematic study of invasive flows as a function of surface tension. Such 
a study would be possible within our model, again by varying the inverse temperature-like parameter p. 
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